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Abstract 

KPW  is  a  time  scale  algorithm  that  combines  Kalman  filtering  with  the  basic 
time  scale  equation  (BTSE).  A  single  Kalman  filter  that  estimates  all  clocks  si¬ 
multaneously  is  used  to  generate  the  BTSE  frequency  estimates,  while  the  BTSE 
weights  are  inversely  proportional  to  the  white  FM  variances  of  the  clocks.  Results 
from  simulated  clock  ensembles  are  compared  to  previous  simulation  results  from 
other  algorithms. 


INTRODUCTION 

The  purpose  of  a  time  scale  is  to  create  a  virtual  clock  from  an  ensemble  of  physical  clocks  whose 
differences  from  each  other  are  measured  at  a  sequence  of  dates  (a  date  being  the  displayed  time 
of  a  clock).  The  virtual  clock  is  defined  as  an  offset  from  one  of  the  clocks,  computed  from  the 
measurement  data  by  some  algorithm.  We  usually  want  the  virtual  clock  to  be  quieter  than  any  of 
the  real  clocks  in  both  the  short  term  and  the  long  term. 

One  approach,  which  was  tried  in  the  early  1980s  [1],  is  to  run  a  Kalman  filter  on  the  clock 
difference  measurements,  the  noise  of  each  clock  having  previously  been  modeled  by  a  stochastic 
linear  system.  The  filter  produces  an  estimate,  unbiased  and  with  minimum  error  variance,  of  the 
phase  and  frequency  of  each  clock;  moreover,  if  we  offset  the  tick  of  each  clock  by  its  phase  estimate, 
we  arrive  at  a  single  point  on  the  time  axis  (if  the  measurements  are  noiseless).  It  makes  sense, 
then,  to  regard  this  point  as  the  estimated  origin  of  the  ensemble,  and  to  use  the  sequence  of  these 
values  as  a  time  scale.  This  time  scale,  which  was  realized  as  TA  (NIST),  was  reported  to  follow 
the  clock  with  the  best  long-term  stability,  regardless  of  its  short-term  stability  [2] .  My  goals  have 
been  to  reproduce  this  finding,  understand  it,  and  improve  the  method.  I  seem  to  have  achieved 
the  first  and  third  goals,  but  not  the  second. 

It  turns  out  that  a  good  time  scale  algorithm  can  be  constructed  by  injecting  some  of  the  Kalman- 
filter  information  into  the  traditional  “basic  time  scale  equation”  (BTSE),  which  requires  frequency 

*This  work  was  performed  by  the  Jet  Propulsion  Laboratory,  California  Institute  of  Technology,  under  a  contract 
with  the  National  Aeronautics  and  Space  Administration. 
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estimates  and  a  set  of  weights.  In  brief,  here  is  the  new  time  scale  algorithm,  called  “Kalman  plus 
weights  (KPW).” 

1.  Initialize  the  Kalman  filter  properly,  and  run  it  on  the  clock  models  and  difference  measurements. 

2.  Throw  out  the  Kalman  phase  estimates. 

3.  Use  the  Kalman  frequency  estimates  in  the  BTSE,  whose  weights  are  made  inversely  proportional 
to  the  white  FM  variances  of  the  clocks. 

After  describing  the  algorithm  in  detail,  I  show  results  from  simulations  of  clocks  with  independent 
white  FM  and  random-walk  FM  (RWFM)  components.  Included  are  comparisons  with  previous 
simulation  results  for  two  other  time  scales  that  make  heavy  use  of  Kalman  filtering:  the  Barnes- 
Allan  “frequency  Kalman”  [3]  and  Stein’s  KAS-21  [4]. 


TERMINOLOGY  AND  NOTATION 

I  shall  try  to  introduce  a  consistent  and  suggestive  notation  in  which  to  work.  The  ensemble  has  n 
clocks  H\y . . .  Hn .  A  date  is  the  displayed  time  of  a  clock,  determined  by  counting  its  oscillations. 

At  date  £,  the  following  quantities  are  defined: 

hi  (t)  =  time  coordinate  (on  some  time  scale)  of  Hl‘> s  tick  when  it  shows  date  t;  not  directly 
observable. 

ho  ( t )  =  time  coordinate  of  an  ideal  clock  Ho;  ho  {t)  —  a  +  bt  for  some  constants  a,  ft.  This  is  not  an 
unattainable  concept;  an  ideal  clock  can  and  will  be  defined  by  extrapolating  the  initial  state  of  a 
physical  clock;  see  the  section  on  startup. 

he  ( t )  =  time  scale  or  ensemble  time,  the  time  coordinate  of  a  virtual  clock  He ,  to  be  determined 
by  its  computed  offsets  from  the  physical  clocks. 

xe  ( t )  =  he  ( t )  —  ho  (£),  offset  of  He  from  Ho ,  also  called  a  time  scale  here. 

Xi  (t)  =  ht  ( t )  —  ho  (£),  offset  of  Ht  from  Ho- 

Xij(t)  =  x%  (t)  —  xj(t)  =  ht  (t)  —  clock  difference  measurements,  taken  at  an  increasing 

sequence  of  dates  to,  ti, _ This  study  assumes  noiseless  measurements. 

Xie  (t)  =  Xi  ( t )  —  xe  ( t )  =  hx  (t)  —  he  (£),  offset  of  Ht  from  He .  The  xie  ( t )  are  to  be  computed  as 
statistics  of  the  measurements  through  date  t,  perhaps  with  some  initial  conditions.  A  time  scale 
xe  ( t )  is  determined  by  Xi  ( t )  and  Xie  ( t )  for  some  i.  When  the  measurements  are  noiseless,  it  usually 
turns  out  that  any  i  can  be  used,  that  is,  xt  ( t )  —  xle  ( t )  gives  the  same  value  xe  (t)  for  all  i.  An 
equivalent  condition  is 

xie  ( t )  -  x3e  ( t )  =  Xij  (t),  *,  j  =  1,  . . .  ,n.  (1) 

If  (1)  is  fulfilled,  I  shall  say  that  the  offsets  xie  (/;)  are  consistent  with  the  measurements.  This  just 

means  that  the  set  of  points  {xie  (t)  :  i  =  1,  ...  ,n}  is  a  rigid  translation  of  the  set  {xi  (t)  :  i  —  1,  . . .  ,n}. 

1 A  tradename  of  Timing  Solutions  Corporation 
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AVERAGE  TIME  SCALES 


This  category  includes  many  of  the  time  scales  in  actual  use  [5].  An  average  time  scale  is  defined 
recursively  at  measurement  date  t  from  measurement  date  t  —  r  by  an  equation  of  form 


n 

xe  ( t )  =  Xe(t-T)  +  'Y^/Wt  ( t )  [Xi  (t)  -xt(t-T)~  T$i  ( t  -  r)] .  (2) 

2=1 

The  weights  Wi  ( t )  (which  add  to  1)  and  the  frequency  estimates  y%  ( t )  depend  on  the  measurements 
through  date  t .  From  (2)  and  the  previous  definitions,  there  follows  the  basic  time  scale  equation 
(BTSE), 


n 

xje(t)  =  22wi(t)[Xji(t)  +  Xie(t-T)  +  Tyt(t-T)\,  j  =  1,  .  . .  ,U,  (3) 

2=1 

a  recursive  equation  for  the  offsets  xje  ( t ).  This  is  the  computation  that  is  actually  performed  to 
obtain  the  offsets  of  the  virtual  clock  from  the  physical  clocks.  These  offsets  are  consistent  with 
the  measurements. 

Average  time  scales  usually  calculate  yt  ( t )  as  an  estimate  of  the  frequency  of  Ht  relative  to  the  scale 
He  as  calculated  through  date  t.  In  the  present  formulation,  yi  (t)  is  an  estimate  of  the  frequency 
of  Hi  relative  to  ffo,  not  He .  This  is  the  case  for  the  Kalman-based  estimate  discussed  below;  the 
Kalman  filter  knows  nothing  about  He . 


CLOCK  MODEL  AND  KALMAN  FILTER 

At  this  stage  of  development,  I  am  using  a  two-state  clock  model:  white  FM  plus  RWFM.  Jones 
and  Tryon  [1]  showed  how  to  integrate  the  differential-equation  model  to  a  stochastic  difference- 
equation  model  for  discrete  measurement  dates,  which  may  be  unequally  spaced.  For  the  ith  clock, 
the  equations  taking  the  state  [#*,  yt\  from  date  t  —  r  to  date  t  can  be  stated  as 


Xi  ( t )  =  Xi(t-r)- f  ry%  (t  -  r)  +  wxi  ( t ,  r) 
Vi  (t)  =  yl(t-r)+  wyi  (t,  r) . 


The  process  noise  vector  [wxl  (£,  r) ,  wyi  (t,  r)]  (uncorrelated  over  dates  and  clocks)  has  covariance 
matrix 


T  0 

t3/3  t2/ 2  ' 

qxi 

0  0 

+  Qyi 

_  t2/2  t 

(5) 
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in  which  the  q’s,  which  we  assume  to  be  known,  specify  the  white  FM  and  RWFM  noise  levels.  In 
terms  of  them,  the  Allan  variance  of  the  ith  clock  is 


(lyiT 

3  ' 


(6) 


The  overall  state  vector  is  X  ( t )  =  [x\  (t)  ,yi(t) .  ,xn  (t) ,  yn  (<)].  In  a  standard  way  [6],  which 
will  not  be  repeated  here,  the  Kalman  filter  uses  the  model  (4)-(5)  and  the  measurements  Xij  ( t )  to 
obtain  a  recursive  estimate  X  ( t )  =  [x\  (t) ,  fa  (t) , . . .  ,xn  (t) ,  yn  (t)]  and  its  error  covariance  matrix 
P  ( t )  from  the  same  quantities  at  date  t  —  r. 

It  turns  out  that  the  Kalman  phase  estimates  X{  (t)  are  consistent  with  the  clock  measurements  in 
the  sense  that  xt  ( t )  —  x3  (t)  =  x^j  (t);  consequently,  it  makes  sense  to  define  a  natural  Kalman  time 
scale  by 


xeK  (t)  =  xt  ( t )  -  Xi  ( t ) 


(7) 


(the  same  for  all  i).  It  is  this  scale  that  was  used  for  TA  (NIST)  and  found  wanting. 

Startup 

The  Kalman  filter  must  be  initialized  at  a  starting  date  t\  by  providing  a  state  estimate  X  (t\) 

By  taking  care  with  this  task,  we  can  establish  a  reference  for  the  ensemble  and  make  the  filter 
settle  down  quickly  [7].  Let  us  take  noiseless  clock  difference  measurements  (to)  ,xtJ  (ti),  where 
t\  =  <o  +  r.  Without  loss  of  generality,  we  can  assume  that  the  X{  (to)  are  known  exactly.  Some 
initial  information  about  the  random  walk  frequency  states  is  needed,  too.  For  this  purpose,  let 
us  regard  Hi  as  a  master  clock  whose  initial  frequency  state,  relative  to  some  ideal  clock  Ho,  can 
be  defined  or  estimated.  Thus,  let  y\  (to)  be  some  unbiased  prior  estimate  of  y\  (to),  with  error 
variance  p\.  We  can  always  set: 


x\  (to)  =  0,  y\  (t0)  =  0,  pi  =  0.  (8) 

The  implication  of  (8)  is  that  Ho  is  defined  as  the  noiseless  extrapolation  of  the  initial  state  of 
Hi,  which,  though  unknown,  can  still  act  as  a  reference.  This  convention  will  be  called  the  Master 
Clock  1  Startup.  This  is  not  the  same  thing  as  using  Hi  as  a  master  clock  during  the  run  [8] ;  after 
startup,  its  state  is  estimated  (relative  to  Ho)  on  a  par  with  the  states  of  the  other  clocks. 

Here  axe  the  initialization  equations,  whose  derivation  is  omitted: 


Vi  (to)  =  2/1  (to)  +  -  [*»i  (h)  -  xa  (t0)] ,  i  =  1,  •  •  ■  ,  n 

T 

Xi  (ti)  =  xt  (to)  +  ryt  (to) ,  yt  (h)  =  fa  (to) ,  i  =  l,...,n. 


These  give  X  (to)  and  X  (ti).  For  the  Master  Clock  1  Startup  (assumed  for  simplicity),  we  have 
P(h)  =  AQ(t)At,  where  Q  (r)  is  the  overall  process  covariance  matrix  with  diagonal  blocks 


448 


Qi  (r),  and 


‘  1  0 
0  1 
1  0  0 


1  0  0 

_  7  0  0 

KPW  ALGORITHM 

The  KPW  time  scale  algorithm  can  now  be  quickly  given.  Set  Xje  (to)  =  X{  {to)  (that  is,  xe  (to)  = 
0).  Initialize  the  Kalman  filter  from  the  clock  difference  measurements  at  to  and  t\.  Run  the 
Kalman  filter  on  the  measurements  at  dates  t2,  t3, . . .  to  produce  state  estimates  (tfc) ,  yt  (tfe).  At 
measurement  date  t  =  t^,  k  >  1,  apply  the  BTSE  (3),  where  t  —  r  =  t-k-i  ■  The  Kalman  phase 
estimates  xt  are  not  used  at  all. 

It  remains  to  specify  the  BTSE  weights.  To  argue  towards  a  reasonable  choice,  define  the  estimated 
RWFM  component  of  Xi  (t)  as  the  discrete  time  integral  of  iji  ( t ), 


0 

1 

0  •••  00 


®I,RWF  (t)  =  Arwf  {t  -  t)  +  ryt{t  —  t) 

^i,RWF  (^o)  =  xi  (to)  . 


(9) 


Assume  constant  weights  wt.  Summing  (2)  over  the  measurement  dates  from  t\  to  t,  we  obtain 


n 

Xe  ( t )  =^2m  [xt  ( t )  -  X2iRWF  (^)]  •  (10) 

1=1 

We  can  regard  the  quantity  in  brackets  as  an  estimate  of  the  white  FM  component  of  Hz.  Thus, 
the  implicitly  defined  time  scale  is  a  weighted  average  of  approximate  white  FM  components.  To 
try  to  minimize  the  instability  of  xe  ( t ),  we  make  w%  proportional  to  l/qXt .  This  method  has  been 
used  for  all  the  simulations.  In  practice,  if  the  q1  s  are  revised  in  the  middle  of  a  run,  then  so  would 
the  w: s. 

SIMULATIONS 

These  were  run  for  various  sets  of  clocks  as  determined  by  their  s.  For  convenience,  the  data 
were  generated  at  equally  spaced  dates.  The  model  equations  (4)  were  used  to  generate  the  true 
clocks.  The  Kalman  filters  were  initialized  by  the  Master  Clock  1  Startup  condition  (8),  and  were 
mechanized  by  a  covariance  square  root  method  [9,10]  to  avoid  numerical  instability  and  problems 
with  singular  covariance  matrices.  Because  the  true  clocks  were  available,  the  time  scales  were 
computed  from  (9)  and  (10)  instead  of  the  BTSE  (3). 
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Two  Opposite  Clocks 

The  simplest  and  most  revealing  example  is  two  clocks  that  are  as  opposite  as  they  can  be:  H\ 
is  pure  white  FM,  while  is  pure  RWFM.  The  results  are  shown  in  Fig.  1  (arbitrary  scaling). 
The  upper  plot  shows  the  phases  of  the  true  clocks  and  the  KPW  time  scale  He.  The  natural 
Kalman  time  scale  (7)  exhibits  an  extreme  form  of  its  well-known  behavior:  it  is  exactly  equal  to 
:r,\  (t).  This  means  that  the  Kalman  phase  estimate  x\  (t)  is  identically  zero;  all  of  His  white  FM 
is  thrown  onto  X2  ( t ). 

The  middle  plot  shows  the  total  frequency  (difference  quotient  of  the  phase)  for  the  KPW  scale 
He  and  the  true  clocks;  the  curves  are  offset  for  clarity.  For  H2,  the  total  frequency  is  the  same  as 
the  RWFM  state  yz  (t).  The  two  upper  intertwined  curves,  which  are  y?  (t)  and  its  estimate  ij2  (t), 
show  that  the  Kalman  filter  does  an  excellent  job  of  estimating  this  frequency  state. 

In  the  lower  plot,  the  straight  lines  are  the  theoretical  Allan  deviations  of  the  true  clocks;  the 
dots  are  the  measured  Allan  deviations.  The  KPW  scale  He  is  only  moderately  noiser  than  H2 
for  short  r,  and  about  the  same  as  Hi  for  long  r.  All  the  weight  is  on  H2;  this  means  that 
xe  (t)  =  X2  (t)  —  ®2,rwf  (£)>  with  H\  playing  no  role  in  the  BTSE.  Nevertheless,  the  long-term 
behavior  of  the  scale  seems  to  be  governed  by  Hi. 

Eleven  NIST  Cs  Clocks 

This  example  comes  from  a  study  by  Barnes  and  Allan  [3]  based  on  simulations  of  a  set  of  cesium 
clocks  whose  q’s  had  previously  been  measured.  Figure  2  shows  the  result  of  a  simulation  using 
the  same  q’s.  The  phases  of  all  the  true  clocks  and  the  total  frequencies  of  three  of  them  are 
shown.  The  crosses  in  the  lower  plot  (data  from  Fig.  9  of  [3]  show  the  Allan  deviation  of  a  time 
scale  derived  from  the  “frequency  Kalman”  filter,  which  uses  pseudo-measurements  of  frequency 
differences.  That  scale  dips  a  little  lower  than  the  KPW  scale  at  r  =  105  s.  At  r  =  104  s,  though, 
the  measured  KPW  Allan  deviation  is  1.30  x  10-14,  close  to  the  minimum  value  1.27  x  10-14  that 
can  be  achieved  by  a  weighted  average  of  the  white  FM  components  of  these  clocks. 


Eight  Imaginary  Clocks 

This  example  reproduces  a  simulation  that  was  carried  out  by  Stein  [4]  on  an  imaginary  eight-clock 
ensemble  to  demonstrate  the  KAS-2  time  scale  algorithm.  The  odd-numbered  clocks  all  have  the 
same  q’s,  as  do  the  even-numbered  clocks.  Figure  3  shows  the  results  of  simulating  this  ensemble; 
the  crosses  in  the  lower  plot  show  the  KAS-2  stability  from  Stein’s  Fig.  1.  For  each  r,  the  measured 
KPW  Allan  deviation  is  less  than  60%  of  the  theoretical  Allan  deviation  of  the  best  clock  for  that 

T. 


CONCLUSIONS 

The  KPW  time  scale  algorithm  has  been  demonstrated  in  a  simulation  playpen  with  perfect  knowl¬ 
edge  of  the  stochastic  clock  models  and  their  noise  levels.  Under  these  conditions,  KPW  seems  to 
be  competitive  with  other  Kalman-based  time  scale  algorithms.  The  natural  Kalman  time  scale, 
which  does  not  use  the  basic  time  scale  equation,  has  again  been  shown  to  be  noisy  in  the  short 
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term  because  the  Kalman  filter  attributes  the  white  FM  noises  to  the  wrong  clocks.  I  do  not  yet 
understand  why  this  happens. 

A  related  symptom  is  the  unbounded  growth  of  portions  of  the  covariance  matrix.  This  growth 
does  not  harm  the  frequency  estimates,  especially  if  the  Kalman  filter  is  mechanized  by  a  square- 
root  method.  Nevertheless,  as  Weiss  and  Weissert  wrote  [2],  “it  is  suggestive  of  an  undesirable 
situation.”  In  view  of  Brown’s  [8]  work,  it  might  be  possible  to  reduce  this  matrix  transparently. 

The  KPW  algorithm  might  serve  as  the  foundation  of  a  time  scale  that  gives  real-time  results.  Of 
course,  a  practical  time  scale  must  provide  for  clock  insertion  and  removal,  outlier  detection  and 
rejection,  jumps  in  phase  and  frequency,  steering,  adaptive  estimation  of  the  qf  s,  and  so  on.  In 
addition,  the  clock  and  measurement  models  should  be  expanded  to  include  random  run  FM,  white 
PM,  and  measurement  noise. 
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Fig.  1.  Two  opposite  clocks.  Clock  1  =  white  FM,  clock  2  =  RWFM 
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